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Signal transduction pathways control various events in mammalian cells such as growth, 
proliferation, differentiation, apoptosis, or migration in response to environmental stimuli. 
Because of their importance, the activity of signaling pathways is controlled by multiple 
modes of positive and negative feedback regulation. Although negative feedback regulation 
primarily functions to stabilize a system, it also becomes a source of emerging oscillations. 
For example, the oscillatory behavior of mitogen-activated protein kinase (MAPK) activity 
has been theoretically proposed earlier and experimentally verified recently. However, the 
physiological function of such oscillatory behavior in biological systems remains unclear. 
To understand the functional aspects of this behavior, one should analyze the oscillation 
dynamics from a mathematical point of view. In this study, we applied the phase reduction 
method to two simple, structurally similar phosphorylation-dephosphorylation cycle mod- 
els with negative feedback loops (Models A and B) and a MAPK cascade model, whose 
dynamics all show oscillation. We found that all three models we tested have a Type II 
phase response. In addition, we found that when a pair of each models A and B coupled 
through a weak diffusion interaction, they could synchronize with a zero phase difference. 
A pair of MAPK cascade models also showed synchronous oscillation, however, when a 
time delay was introduced into the coupling, it showed an asynchronous response. These 
results imply that structurally similar or even identical biological oscillators can produce dif- 
ferentiated dynamics in response to external perturbations when the cellular environment 
is altered. Synchronous or asynchronous oscillation may add strength to or dampen the 
efficiency of signal propagation, depending on subcellular distances and cell density. Phase 
response analysis allows prediction of dynamics changes in oscillations in multiple cellular 
environments. 
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1. INTRODUCTION 

Oscillatory dynamics are widely distributed in nature (Strogatz, 
2003) . In biological systems, circadian rhythm, heart rhythm, loco- 
motion, and electrical activity of the brain are well-known oscil- 
lation generators (Winfree, 2001; Khalsa et al, 2003). Although 
oscillatory behavior is a product of negative feedback regulation, 
the question of how the oscillatory information is processed in 
biological systems is still unresolved. Mammalian cells respond to 
extra- cellular signals and transfer this information to the nucleus 
to express/repress genes necessary for adaptation to a new envi- 
ronment or differentiation state. Signal transduction pathways 
play important roles to control expression of the correct genes 
and with the precise timing to satisfy cellular needs. Therefore, 
signaling pathways are spatio-temporally controlled by many pos- 
itive and negative feedback loops through transcriptional and 
post-transcriptional modification. As a result, several types of 
oscillatory behaviors in the components of signaling pathways 
can be observed when negative feedback regulation is introduced 
into a system. For example, Shankaran et al. have shown persis- 
tent periodic shuffling of fluorescent-labeled extra-cellular signal- 
regulated kinase (ERK, a subset of the family of mitogen-activated 



protein kinases, MAPK) between cytosol and nucleus in epider- 
mal growth factor (EGF) stimulated cells at the single cell level. 
Intriguingly, these periodic cycles among neighboring cells were 
asynchronous (Shankaran et al., 2009). ERK is one of the deter- 
ministic kinases that control transcription when translocated into 
the nucleus, therefore this nuclear shuffling process is highly reg- 
ulated. The work of Shankaran et al. was the first experimental 
demonstration of the oscillatory behavior of ERK, although this 
was predicted earlier on theoretical grounds (Kholodenko, 2006). 
Given this asynchronous oscillation, one would think that it would 
be difficult to identify ERK dynamics in a population of cells, where 
the signal would be averaged. In addition, since periodic activation 
of ERK has been difficult to demonstrate experimentally, the cellu- 
lar conditions leading to oscillatory ERK activation are likely quite 
narrow and restricted. In general, when oscillators interact with 
each other through a strong coupling, they tend to synchronize. 
Therefore, the asynchronous oscillation observed by Shankaran 
et al. suggests that the coupling strength of ERK in neighbor- 
ing cells is weak, at least under the experimental conditions used 
in these studies. The question then arises, what kind of condi- 
tions allows a pair of cells to achieve asynchronous oscillation? Is 
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the weak coupling enough to cause asynchronous oscillation? We 
have investigated these questions by applying the phase reduction 
method to two models of phosphorylation-dephosphorylation 
cycles and in a MAPK cascade, all of which exhibit negative 
feedback regulation. 

The mechanism that causes the emergence of oscillations 
has been energetically studied (for example, Guckenheimer and 
Holmes, 1983; Kuramoto, 2003). Theoretically, oscillators are clas- 
sified according to their bifurcation types, such as saddle-node 
and Hopf bifurcations. Oscillation by the saddle-node bifurca- 
tion emerges when a half-stable cycle splits into a pair of stable 
and unstable limit cycles, but oscillation by the Hopf bifurcation 
emerges when a stable spiral fixed point changes to a unstable 
spiral fixed point surrounded by a stable limit cycle (Hale and 
Kocak, 1991; Strogatz, 2001). Although phase space structures can 
be partly derived from such bifurcation types, the dynamic prop- 
erties of the system have to be evaluated by other methods. To 
investigate the underlying oscillatory mechanism, a framework 
termed the phase reduction method has been developed in math- 
ematics and non-linear physics (Hansel et al., 1995; Hoppensteadt 
and Izhikevich, 1997; Kuramoto, 2003). By using this method, an 
oscillator in a high dimensional space can be described by only 
one variable, phase, and its dynamics are packed into a phase 
response (or phase sensitivity) function. The phase response func- 
tion has been derived analytically, not only from mathematical 
models but also from experimental biological data (Reyes and Fetz, 
1993; Khalsa et al, 2003; Lahav et al, 2004; Strieker et al, 2008). 
This method facilitates the classification of structurally related 
but dynamically differentiated biochemical oscillators. The phase 
response functions have been classified into two types, commonly 
referred to as Type I and Type II (Hansel et al., 1995; Hoppensteadt 
and Izhikevich, 1997; Kuramoto, 2003). A Type I phase response 
function generally attains a positive value during an oscillation 
period, whereas a Type II phase response function possesses a sig- 
nificantly large region of negative values. A small perturbation to 
an oscillator advances its phase when it is in a phase that gen- 
erates a positive phase response, but retards its phase when in 
a negative phase response. It is known that Type I and Type II 
phase response functions correspond to saddle-node and Hopf 
bifurcations, respectively (Rinzel and Ermentrout, 1998). 

In this study we have used three models and the phase reduc- 
tion method to investigate the type of phase response function and 
phase difference of two weakly coupled oscillators in the steady 
state. 

2. RESULTS 

2.1. SINGLE OSCILLATOR: A PH0SPH0RYLATI0N- 
DEPH0SPH0RYLATI0N CYCLE 

Several modes of negative feedback regulation have been identified 
in signal transduction pathways, and these are potential candidates 
for emerging oscillatory phenomena. Many fundamental negative 
feedback models that cause the emergence of oscillations have been 
proposed (Kholodenko, 2006; Novak and Tyson, 2008). Here, we 
adopt the simple phosphorylation-dephosphorylation cycle mod- 
els proposed by Kholodenko (2006). While he has considered all 
the possible topologies of feedback regulation to phosphorylation 
and dephosphorylation steps in the cycle, we use two, Models A 



and B (the latter of which corresponds to Model C in the original 
paper) in our study (Figures 1A,F, Section 4.2). In these models, 
negative feedback is realized by inhibiting kinase (Kin) production 
(or its activity) in Model A and enhancing phosphatase (Phos) 
production (or its activity) in Model B. First, we explore the para- 
meters exhibiting oscillation by varying Phos and Kin for Models 
A and B, respectively (Figure 1). The parameter regions that can 
induce oscillations, resultant oscillation periods, and frequencies 
are shown in Figures 1B,G. The long-dashed lines indicate the 
parameter values that we have adopted in the following analyses. 
Figures 1C,H show the periodic orbits of the models. The oscilla- 
tion periods in the two models are clearly very different from each 
other. 

We next applied a phase reduction method to these models 
and calculated the phase response functions of Models A and B. 
The details of phase reduction method are found in Section 4.1, 
in which the state vector X(t) is given by (Mp(f), Kin(f)) and 
(Mp( t), Phos( f ) ), where Mp represents an activated and phospho- 
rylated form of M, for Models A and B, respectively. Figures 1D,I 
show the periodic orbits for one oscillation period in Models A 
and B, respectively, in which the peaks of the Mp (blue lines) 
are located at the origin of the phase. The phase response func- 
tions of Mp of both models are similar to each other and have 
significantly large regions of both positive and negative values 
(Figures IE, J for Models A and B, respectively), which means that 
they can be categorized as Type II oscillators. This result implies 
that these regulatory networks have similar phase responses of 
Mp to a small perturbation regardless of the difference in their 
biological feedback targets. 

2.2. COUPLED OSCILLATORS: INTERACTING 

PH0SPH0RYLATI0N-DEPH0SPH0RYLATI0N CYCLES 

Next, we investigated the behavior of the above models in the 
presence of a weak interaction in the steady state by presum- 
ing that two cells are located next to each other. The phase 
reduction method allows calculation of the fixed phase differ- 
ence of a weakly coupled pair of identical oscillators in the 
steady state (for details, see Section 4.1). Here, we consider the 
case of two identical oscillators interacting through a weak dif- 
fusion coupling (see Figures 2A,C). In addition to the model 
equations for Models A and B (Section 4.2), we adopt the fol- 
lowing diffusion couplings. For Model A, the interaction func- 
tion in Equations (7) and (8) is given by C((Mp,-(f), Kin,(t)), 
(M Vj (t),nnj(t))) = (-gMp(Mp,(f) - Mp;(f)), -sKin(Kin,(f) - 
Kinj(f))). Similarly, for Model B, it is given by C((Mp,(f), 
Phosf(r)), (M P; (f ), Pho Sj -(r ))) = (-^Mp(Mp,(f ) - M Pi (t)), 
— gPhos(Phos,(f) — Phosj(t))). We assume that the coefficients of 
interactions, gMp, gKin, and gPhos, are sufficiently small for each 
oscillator to remain in the basin of the periodic orbit. Under weak 
coupling conditions, a coupled dynamical system can be reduced 
to a system governed only by phase difference. The reduced 
dynamical system is given by Equations (17) and (18). Once we 
obtain the phase response function, we can calculate the gamma 
function (T~((/>), defined by Equations (17) and (18)), thereby 
delineating the dynamics of the phase difference. The gamma 
functions for Models A and B are shown in Figures 2B,D, respec- 
tively. As shown in the figures, (/> = 0 and 0.5 satisfy r _ (0) = 0 
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FIGURE 1 | Phosphorylation-dephosphorylation cycle models. (A) same as in (C). (E) Phase response functions of the model, which are 

Schematic of Model A, in which the phosphorylated M, Mp, decreases the obtained by solving the adjoint equation given in Equation 4 (Ermentrout, 

production of kinase or inhibits its activity. Here, Phos takes a constant value, 1996). Colors are the same as in (C). (F) Schematic of Model B, in which the 

200. (B) Period (blue line) and frequency (green line) of oscillation as a phosphorylated M, Mp, increases the production of phosphatase or activates 

function of Phos. (C) Oscillatory behavior of the model network. Blue and red its activity. Here, Kin takes a constant value, 1 50. (G-J) Same as in (B-E), but 

lines correspond to Mp and Kin, respectively. (D) Periodic orbit during one here they are obtained using model B. Blue and red lines in (H-J) correspond 

period of the model. The horizontal axis represents the phase. Colors are the to Mp and Phos, respectively. 



for both models, and only <p = Q satisfies the stability conditions 
(Equations (19) and (20)). Therefore, the oscillators in both cou- 
pled dynamical systems are asymptotically synchronized in the 
steady state. Next, we evaluated the shape of the gamma function 
by varying the ratio gKin/gMp (gPhos/gMp) between 0 and 1 for 
Model A (Model B) and find that 4> = 0 is the only stable solution 
in this parameter region. Therefore, these results suggest that in 
a many body system consisting of each model, their oscillation 
cycles can synchronize with an almost zero phase difference in a 
noisy environment. 

2.3. MAPK CASCADE MODEL 

Next, we considered the asynchronous oscillations experimentally 
observed by Shankaran et al. (2009). They reported that EGF 
induces oscillations in the nuclear localization (an indication of 
activation) of ERK in living cells and that the oscillations are asyn- 
chronous between neighboring cells (Figure 1 in Shankaran et al, 
2009). Here, we adopt the MAPK cascade model that was origi- 
nally developed by Huang and Ferrell (1996) and later modified 
by Qiao et al. (2007). The Huang-Ferrell model has been widely 
used for analyzing the dynamic behavior of the MAPK system 
(for example, Ferrell and Machleder, 1998). The original model 
and parameters produce an ultrasensitive MAPK activity with a 
Hill coefficient of 4-5, but later, Qiao et al. found that the model 
can also produce oscillations with other parameter sets in a wide 



range of the parameter space. Here, we used the Huang-Ferrell 
model modified by Qiao et al. to analyze the oscillatory behavior 
of MAPK activity (Figure 3A, Section 4.3). 

The time evolutions of MAPKKPP and MAPKPP in an oscilla- 
tory state are shown in Figure 3B. Here, we assume that MAPKPP 
is an activated ERK as observed by Shankaran et al. We applied the 
phase reduction method to this model and evaluated the stable 
fixed phase difference by calculating the phase response function 
of the model. The phase response functions of MAPKKPP and 
MAPKPP are shown in Figures 3C,D, respectively, in which the 
origin of the phase corresponds to the peak of MAPKPP. As shown 
in the Figure 3D, MAPKPP showed a Type II phase response, the 
same as those of Models A and B described in the previous section. 

Next, we considered the phase difference of a pair of iden- 
tical MAPK oscillators. Because the entire signal transduction 
network in the cells used by Shankaran et al. remains unknown, we 
substituted a MAPK cascade as an EGF-induced signaling path- 
way in the observed neighboring cells. In addition, we assumed 
that the interaction between a pair of cells can be effectively 
modeled as a function of the difference between the MAPKPPs. 
Then we adopted a simple diffusion type interaction function: - 
gMAPKPP(MAPKPP,(t) - MAPKPPyU))- Using this function, we 
calculated the gamma function as shown in Figure 3E (black line). 
The value <p = 0 is the only stable solution and the two MAPK 
oscillators are synchronized asymptotically. Thus, employment of 
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FIGURE 2 | (A) Schematic of a pair of the identical oscillators of Model A. 
Here, we assume weak interactions between each Mp and Kin. (B) Gamma 
function of the coupled system shown in (A), in which we adopt gMp = 1 and 



gKin = 1. (C) Same as (A), but Model B was used by assuming weak 
interactions between each Mp and Kin. (D) Gamma function of the coupled 
system shown in (C), in which we adopted gMp = 1 and gPhos= 1. 




FIGURE 3 | (A) Schematic of the Huang-Ferrell model. Here, 
MAPKKK* indicates the activation form of MAPKKK. (B) Oscillatory 
behaviors of MAPKPP and MAPKPP in the steady state of the 
model. (C,D) Phase response functions of MAPKPP (C) and 



MAPKPP (D), respectively. (E) Gamma function without a time 
delay (black line) and one with a time delay t = 0.25 (red line). (F) 
Stable (blue) and unstable (red) fixed phase differences as a function 
of time delay r. 
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only a diffusion interaction between the two MAPK models (or 
cells) was insufficient to reproduce the asynchronous oscillation 
of MAPKPP that has been detected in living cells. 

To achieve an asynchronous oscillation, the sign of the coef- 
ficient can be changed. When doing so, we obtained a gamma 
function whose shape is the reflection of the one shown as a 
black line in Figure 3E and <f> = 0.5 results in a stable solution. 
Because we considered an effective interaction in our model, a 
negative diffusion coefficient can be considered. Another possible 
scenario is to incorporate a time delay within the interaction, - 
gMAPKPP(MAPKPP,-(f) - MAPKPP ; (t - t)), where r > 0. This 
scenario is reasonable because spatially isolated cells need to com- 
municate with each other, but the signal from one cell will be 
delivered with a time delay to another cell when the two cells 
remain apart. In agreement with this theoretical assumption, oscil- 
latory ERK activity was only observed in cells grown at low density 
but not in cells at confluency (Shankaran et al., 2009). When we 
incorporated a time delay into the interaction, the gamma func- 
tion could be very simply calculated (Equations (30) and (31)). 
Considering x = 0.25 as a delay time, we calculated the gamma 
function in Equation (31), as illustrated in Figure 3E (red line). As 
a result, <p = 0.5 transforms into a stable solution. Therefore, ERK 
signals potentially asynchronously oscillate in the delayed system. 
In Figure 3F, we show how stable and unstable fixed phase solu- 
tions alter by varying the time delay, r, between 0 and 1. A time 
delay 0.24 < r < 0.63 results in <p = 0.5, which is the only stable 
fixed phase difference. Interestingly, a wide range of time delays 
could induce asynchronous oscillations of MAPK. 

3. DISCUSSION 

In this study, we first evaluated the phase response function, 
a characteristic property related to oscillation, for two simple 
phosphorylation-dephosphorylation cycle models A and B, in 
which negative feedback regulation either inhibits the kinase 
activity or enhances the phosphatase activity. The two models 
have relatively different oscillatory periods. However, their phase 
response functions corresponding to the reaction product Mp are 
very similar and both have a Type II phase response, i.e., their 
phase can be retarded or advanced depending on the timing of 
an external stimulus. This result suggests that even if cells use 
different negative feedback regulatory mechanisms (e.g., kinase 
inhibition or phosphatase activation), they can produce simi- 
lar, although not exactly the same, oscillatory dynamics. It also 
indicates that similar oscillatory signals can be generated using dif- 
ferent biological components and resources, such as kinases and 
phosphatases, produced by regulated transcription or activated by 
post-transcriptional modification in a cell context manner. Fur- 
thermore, we investigated the fixed phase differences of a pair 
of identical oscillators interacting through weak diffusion cou- 
pling. Here we found that the oscillations originating from each 
model can synchronize even if their interaction is weak. This result 
implies that when two almost identical oscillations (i.e., the bio- 
chemical reaction dynamics generated from Model A or B) can 
interact in spatially separated subcellular locations (such as the 
cytoplasm and nucleus), their oscillatory signals can still be syn- 
chronized. These results suggest that oscillation dynamics can be 
robustly maintained within the cell. 



We next analyzed the MAPK cascade and found that the phase 
response function corresponding to MAPKPP is also Type II. 
When two identical MAPK oscillators are coupled through dif- 
fusion, these signals are synchronized. Again, in addition to the 
above two models, this result confirmed that oscillatory signals can 
be robust within the cell. However, the presumed interaction of 
MAPKs mediated by diffusion in our analysis did not satisfactorily 
reproduce the asynchronous ERK oscillation dynamics reported by 
Shankaran et al. (2009). We reasoned that this inconsistency could 
be the result of a time delay in the coupling and, under this condi- 
tion, asynchronous oscillation was observed with a wide range of 
parameters in time delays. Our results suggest that whether a pair 
of cells oscillates synchronously or asynchronously may depend 
on the distance between the cells. In the experimental setting of 
Shankaran et al., oscillatory ERK activity was only observed in 
cells grown at low density but not in confluent cells. Therefore, we 
presume that compensatory mechanisms would mask the oscilla- 
tory dynamics of ERK in cells at high density. In reality, thousands 
of MAPK molecules exist in a cell, therefore molecular regula- 
tion of MAPK in a cell in a tightly packed cell population may be 
quite different from that in a sparse cell population, and there- 
fore the former circumstance would interfere with experimental 
visualization of oscillatory dynamics of MAPK. 

In this study, we only dealt with identical pathways under weak 
interaction to examine the effect of two units coupling for an oscil- 
latory response. However, our current study suggests that a pair 
of slightly different oscillators under weak interaction conditions 
can result in similar synchronous behavior: therefore the method 
should be applicable for an evaluation of non-identical oscillatory 
units as well. The phase reduction method thus can be applied 
to pairs of widely different oscillator modules, as are frequently 
observed in biological networks. 

4. MATERIALS AND METHODS 
4.1. PHASE REDUCTION METHOD 

The phase reduction method has been widely applied to coupled 
oscillator systems. Here, we briefly describe the derivation of the 
phase model. First, we define the phase. We consider a system 
whose dynamics is described in the vector form as follows: 



dX(t) 
dt 



F(X(t)), 



(1) 



where X(t) represents the dynamical variable, and F(X(t)) is a 
vector function that determines the dynamics of the system. Here- 
after, we restrict the phase space of the system to a region in its 
basin of attraction in which a stable limit cycle, /(f), exists. We 
assume that the period of the solution is T, i.e., x (t+ T) = x(t). 
According to the phase reduction method, such a system can be 
reduced to a system consisting of the phase degree of freedom only. 
We can define a phase <p{t) e [0, 1) along the periodic orbit %(t) 
with constant time derivative as follows: 



d<t>(t) 
dt 



1, 



(2) 



in which we can define the origin arbitrarily. Because the periodic 
orbit can be parameterized by <p{t), we describe it by x(</>( f ))- 
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Although the phase is defined along the periodic orbit, it can be 
extended into the neighborhood of the periodic orbit, which is 
known as the asymptotic phase. Because the asymptotic phase 
can be described as a function of the state X(t), we can derive its 
dynamical equation as follows: 



d<p(X(t)) _ d<p(X) dX(t) 
Jt ~ ~~9~X dT~ 
90 (X) 



(3) 



where • denotes the inner product. From Equations (2) and (3), 
4>(X{t)) must satisfy the following expression 



90 (X) 
dX 



F(X) = 1. 



(4) 



Next, we derive the phase response of the system. We assume 
that the system Equation (1) receives a weak perturbation: 



dX(t) 
dt 



F(X(t)) + p(t). 



(5) 



Similar to Equation (3), the asymptotic phase obeys the 
following dynamical equation: 

d(j>(X(t)) _ 90 (X) dX(t) 
It ~ ~~9X ~dT~ 

-^•<wo> + *o) m 

= i + Z(<p) ■ p(t), 

in which Z(0) is called as the "phase response function" or "phase 
sensitivity." The phase response function describes the phase 
response (phase shift) of the system to small perturbations. 

Now, we consider a coupled oscillator system consisting of two 
identical oscillatory units, 



dXi(t) 

dt 
dX 2 (t) 

dt 



= F(X l (t)) + C(X l (t),X 2 (t)), (7) 
= F(X 2 (t)) + C(X 2 (t),X 1 (t)). (8) 



Here, X,(f) represents the dynamical variables corresponding 
to the i-th unit, and Cis the interaction function. When the inter- 
action is weak ("weak" means that the units do not leave the basin 
of attraction of the periodic orbit), the coupled system can be 
described by the asymptotic phase as follows: 



dMt) 

dt 
dMt) 

dt 



l + Z(0i(O)-C(0 1 (f),0 2 (f)), (9) 

i + Z(<j> 2 (t)) -C(<pAt), do) 



By replacing the phase 01.2(f) with the phase 01.2(f) defined 
by 0i,2(f) = f+ 01.2(f), Equations (9) and (10) are transformed 
into the following equations: 



dfi(t) 

dt 
djf 2 (t) 

dt 



r(0i(f)- 02(f)), 
r(0 2 (O- 0i(O), 



(ii) 

(12) 



where 



r(0,(f)-0 ; (O)= f 1 ^z(e + 0,( f ))-c(f5 + 0 J (f),e + 0 J (f)), 

Jo 

(13) 

in which we used the fact that 0 1.2(f) vary slowly, the weak 
interaction assumption, and the average during one period. By 
substituting 01.2(f) = — f+0i,2(f) into Equations (11-13), we 
obtain the phase description of the coupled oscillator system as 
follows: 



d<pi(t) 

dt 
dMt) 

dt 



= i + r(0! (t), 02(f)), 
= i + r(0 2 (t), 0i(f)), 



(14) 
(15) 



r(0 



«(f),0;(f))= / dez(0 + <Pi(t))-c(e + <Pi(t),e + <Pj(t)). 

Jo 

(16) 



By subtracting Equation (15) from Equation (14), the ordinary 
differential equation (ODE) is obtained that governs the evolution 
of the phase difference cf>(t) = (f> l (t) — 02(f), 



d(j>(t) 
dt 



= r-(0(O) = r(0(O)-r(-0(O), (17) 



T(0(f))= [ d6Z(6+cP(t))-C(6+<P(t),6). (18) 
Jo 

The fixed phase difference 0 is stable if it satisfies the following 
conditions: 



r-(0) = o, 
dr-(<j>) 



d(j> 



< 0, 



and it is unstable if it satisfies the following conditions: 

r-(0) = o, 
dr-(0) 



rf0 



> o. 



(19) 
(20) 

(21) 
(22) 



Finally, we derive a phase model corresponding to a coupled 
oscillator system with a delayed interaction. In this case, the system 
is expressed as follows: 



dX x {t) 

dt 
dX 2 (t) 

dt 



F(X 1 (f)) + C(X 1 (f),X 2 (f-r)), (23) 
F(X 2 (t)) + C(X 2 {t),X 1 (t-r)). (24) 
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Then, Equations (9) and (10) are changed into the following 
equations: 



The parameters of Models A and B are given below. 



d<pi(t) 
dt 

d<p 2 (t) 
dt 



= 1 + Z(<p 1 (t))-C(<p 1 (t),<j> 2 (t-T)), (25) 
= 1 + Z(0 2 (f))- C(Mt),Mt ~ r)). (26) 



In a first order approximation, we can assume that 
<Pi,i(t — t) = <pi, 2 (t) — x. Then, the phase description of the 
system is given by the following equations: 



d<t>i(t) 

dt 
d<p 2 (t) 

dt 



r(0, 



= 1 + TC0! (0,02(0), (27) 

= i + r(0 2 (o,</>i(O), (28) 
dez(0 + 4>i(t))-c(e + 4>i(t),e 



i{t),4>j{t)) = f 

Jo 

+ (pj(t)-T) 



(29) 



Therefore, the ODE of the phase difference attains the following 
form: 



d<j>(t) 
dt 



= r-(cj>(t)) = r(0(t)) -r(-</>(0), 



(30) 



/' 

Jo 



r(<P(t))= I d8Z(9 + (P(t))-C(9 + (P(t),9-r), (31) 
Jo 

in which the delay is no longer explicitly included in the equation. 

4.2. SIMPLE PHOSPHORYLATION MODEL 

Model A 



dMp 
dt 



dKin 
~dT 



n,n - Vphos 

KZKinM l+AMp/K a ^[PhosMp 
K m \ + M ' 1 + Mp/K a K m2 + Mp ' 

synth deg 



o 1 + Mp/K { _ d e g 
V kin - • - K «~« 



1 + I ■ Mp/K, 



kin 



Kin, 



M = M m - Mp. 

Model B 

dMp _ 

^ — v kin ~ v phos 



k^KinM \+AMp/K a ^PhosMp 
Kmi + M 1 + Mp/K a K m2 + Mp ' 



dPhos 



synth dee 

= V r — V , 

^ phos phos 

_ o i+A dp Mp/K d de 

~ V Phos i + M p/K d V hos > 

M = M m - Mp. 



(32) 

(33) 

(34) 
(35) 
(36) 

(37) 

(38) 

(39) 
(40) 
(41) 



Parameter 


Model A 


Model B 


M tot 


300.0 


300.0 


/.cat 
kin 


1.0 


1.0 


A 


100.0 


100.0 


Ka 


500.0 


500.0 


Km] 


500.0 


500.0 


heat 
*phos 


1.0 


1.0 


Km2 


10.0 


10.0 


Phos(A), Kin(B) 


200.0 


150.0 


^CA). V° h0S (B) 


150.0 


200.0 


K/(A), K d (B) 


100.0 


100.0 


/(A), A dp (B) 


7.5 


7.5 




1.0 


1.0 



4.3. HUANG-FERRELL MODEL 

The ODEs of the model are as follows: 



dv 0 
dt 
dv\ 

dt 

dv 2 

dt 

dv 3 

dt 

dv 4 

dt 
dvs 

dt 
dv 6 

dt 
dv-j 

dt 

dvt, 
dt 

dvg 

dt 
dvio 
dt 

dvn 

dt 
dv\ 2 

dt 
dvn 

dt 
dvu 

dt 



■ oqCqci - (d 0 + ko)v 0 , (42) 

: k>v 0 - a\V\c 2 + d\v 2 - a 2 vic 3 + (k 2 + d 2 )v 3 

- a 4 v 4 vi + (k 4 + d 4 )v 6 , (43) 
: aivic 2 - (d x + fci)v 2 , (44) 

: a 2 c 3 vi - (d 2 + k 2 )v 3 , (45) 

: k 2 v 3 - a 3 v 4 c 4 + d 3 v 5 - a 4 v 4 v x + d 4 v 6 + k 5 v s , (46) 

: a 3 v 4 c 4 - (d 3 + k 3 )v 5 , 

■ a 4 v 4 v\ - (d 4 + k 4 )v 6 , 

■■ k 4 v 6 - a 5 vjc 4 + d 3 v% - a 6 Y]C 3 + (d 6 + /c6)v 9 

- a%v w v 7 + (d & + fc 8 )vi2, 
: a 5 v 7 c 4 - (d 5 + fc 5 )v 8 , 

: aec 5 v 7 -(d 6 + 



(47) 
(48) 

(49) 
(50) 

(51) 



Afev 9 - a 7 v w c 6 + d 7 v u - a & v l0 v 7 + d & v X2 + fc 9 v 14 , 



a 7 v w c 6 - (d 7 + k 7 )v n , 
fl 8 vi 0 v 7 - (d s + k s )v u , 
hv u - a 9 v l3 c 6 + d g v u , 
a 9 vi 3 c 6 - (d 9 + k 9 )vi 4 . 



(52) 
(53) 

(54) 

(55) 

(56) 
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The algebraic equations and total quantities of the molecules 
are as follows: 



CO 
c\ 

C2 
C3 

CA 
C5 
C6 

KKK to t 
El tot 
E2 tot 
KK tot 
KKP'ase tot 

Ktot 
KP'ase tot 



: KKKtot - (VQ + VI + V2 + v 3 + v 6 ), 

: El to t - v 0 , 
■■ E2 tot - v 5 , 

: KK tot - (v 3 + v 4 + v 5 + v 6 + v 7 

+ V 8 + Vg + V12), 

: KKP'asetot - (v 5 + vg), 
: K tot - O9 + v w + Vn + V12 + V13 + v u ), 
■■ KP'ase to t - On + v u ), 
003, 
007, 
004, 



3.7112£ - 
1.0000£ - 
1.2115£ - 
4.0658, 
9.4579£ - 
2.0, 

7.5372£ - 



005, 



002. 



The kinetic parameters are given below: 



(57) 
(58) 
(59) 

(60) 
(61) 
(62) 
(63) 
(64) 
(65) 
(66) 
(67) 
(68) 
(69) 
(70) 



so 


= 1.09625+003 


do 


= 4.32075+001 


ko 


= 6.95255+001 


ai 


= 1.55765+003 


d. 


= 9.43945+001 


*1 


= 2.88745+002 


32 


= 1.91795+003 


d 2 


= 7.52165+001 


k 2 


= 4.34325+001 


33 


= 3.68945+002 


d 3 


= 4.17105+002 


k 3 


= 7.15055+002 


34 


= 4.58365+003 


di 


= 3.37425+002 




= 1.69575+002 


35 


= 2.02195+003 


ds 


= 4.39055+002 


ks 


= 3.48425+002 



a 6 = 2.66345+ 003 


de 


= 5.95985+001 


k 6 = 4.19545+001 


37 = 1.44355+ 003 


di 


= 4.01015+001 


k 7 = 6.34335+ 001 


a 8 = 5.43665+002 


da 


= 2.22115+002 


^8 = 1.17165+002 


39 = 4.3990 5 + 002 


d 3 


= 1.76425+ 002 


k 9 = 4.67055+ 001 


Molecular names, variables, and initial conditions. 


Molecular name 




Variable 


Initial condition 


KKK ■ 5, 




vo 


3.38755-009 


KKK* 




v-i 


1.65615-006 


KKK* ■ E 2 




v 2 


8.15695-010 


KK - KKK* 




v 0 


1.00565-004 


KKP 




v 4 


2.15685-001 


KKP ■ KKP'ase 




V5 


6.10785-006 


KKP ■ KKK* 




V6 


3.22935-006 


KKPP 




V7 


7.04395-003 


KKPP - KKP'ase 




V8 


1.57165-006 


K ■ KKPP 




^9 


8.10415-002 


KP 




V10 


4.94315-001 


KP-KP'ase 




V11 


5.36005-002 


KP-KKPP 




V 12 


5.57955-003 


KPP 




V13 


9.12795-001 


KPP-KP'ase 




V U 


1.39965-002 
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